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Recent progress on studies of the nanoscale mechanical responses in disordered systems has highlighted 
a strong degree of heterogeneity in the elastic moduli. In this contribution, using computer simulations, 
we study the elastic heterogeneities in athermal amorphous solids, composed of isotropic, static, sphere 
packings, near the jamming transition. We employ techniques, based on linear response methods, that are 
amenable to experimentation. We find that the local elastic moduli are randomly distributed in space and 
are described by Gaussian probability distributions, thereby lacking any significant spatial correlations, 
that persists all the way down to the transition point. However, the shear modulus fluctuations grow as 
the jamming threshold is approached, which is characterized by a new power-law scaling. Through this 
diverging behavior we are able to identify a characteristic length scale, associated with shear modulus 
heterogeneities, that distinguishes between bulk and local elastic responses. 

PACS numbers: 83.80.Fg, 61.43.Dq, 62.25.-g 


When traditional, crystalline solids are linearly de¬ 
formed, their elastic responses are typically described by 
affine deformations J!]]. Contrary to this, disordered solids, 
such as thermal amorphous solids, i.e. glasses, disordered 
crystals 10], as well as athermal jammed solids 10], exhibit 
strongly non-affine responses to elastic deformations. This 
non-affine character becomes significantly apparent during 
shear deformation 10] ■ Under shear, constituent particles 
undergo additional non-affine displacemenfs IH], leading 
fo a decrease in fhe shear modulus from a value predicfed 
by fhe affine response only 10] ■ U is Ibis non-affine char- 
acfer fhaf dominafes fhe shear modulus on approach fo fhe 
jammin g fra nsifion, where a mechanically sfable solid loses 

rigidifySBl- 

The appearance of non-affine response is closely relafed 
fo elasfic heferogeneifies i^, especially spafially varying 
shear moduli. Indeed, DiDonna and Lubensky |0] pro¬ 
posed thaf non-affine displacemenfs of parficles subjecf fo 
shearing are driven by randomly flucfuafing local elasfic 
moduli. Amorphous solids refiecl such inhomogeneous be¬ 
havior in fheir mechanical responses af fhe nanoscale [[13- 
1 j], as seen in bofh compufer simulafions fl^ and exper- 


imenfs Gl. Manning and co-workers jislil^ identified 
soff spofs as regions of afypically large displacemenfs in 
low-frequency, quasi-localized vibrational modes. Parfi- 
cle rearrangemenfs, acfivafe d by mechanical load Emi 
and by thermal energy |T^ 18], are therefore understood 
to be spatially correlated with those soft spots, which can 
be linked to locally unstable regions with negative shear 
moduli 0. Furthermore, Ellenbroek et al. l[l3l demon¬ 
strated that the elastic response of jammed packings to lo¬ 
cal forcing fluctuates over a length scale Independently 
Eemer et al. li^ showed that the local elasticity is gov¬ 
erned a different length Recently Karimi and Mal¬ 
oney ll^ reconciled these differing views by considering 
the behaviors of longitudinal and transverse components of 
elastic response. 


Thus, it appears that spatial heterogeneities in local elas¬ 
tic moduli are a key feature to understanding mechanical 
properties of disordered solids. In this contribution, we 
study the elastic heterogeneities in athermal jammed solids 
close to the jamming transition. Specifically, we address 
fhe following poinfs: (i) How are fhe local elasfic mod¬ 
uli distribufed in space? (ii) How do fhose disfribufions 
evolve on approach fo fhe jamming fransifion? (iii) Is fhere 
a lengfh scale over which fhe local elastic moduli flucfuafe? 
For afhermal systems sfudied here, fhe packing fraction cj) 
acis as a confrol parameter fhaf we use fo systematically 
probe sfafic packings of varying rigidify. We characterize 
rigidify by fhe disfance, A(/> = cj) — (pc, from fhe fransi- 
fion poinf pc, or equivalenfly the packing pressure, p. The 
approach of pc from above (p -A O'*") is governed by var¬ 
ious power-law scaliims with Ap in quantities including 
global elastic moduli 0,0, [I^j- In the following, we unveil 
new power-law scalings in the spatial fluctuations of elastic 
moduli. 


Our numerical system consists of N monodisperse, fric¬ 
tionless spheres of diameter a and mass m, in three di¬ 
mensional, periodic, cubic simulation boxes ll23l . Parti¬ 
cles interact via a finite-range, purely repulsive potential; 
V (r) = (e/a)(l — rjaY for r < a, otherwise V (r) = 0, 
where r is the center-to-center separation between two par¬ 
ticles. Here, we show results only for Hertzian contacts, 
a = 2.5 i23n . Eength, mass, and time are presented in 
units of CJ, m, and r = [ma^ . We prepared sys¬ 
tems over several orders of magnitude in packing pressure, 
10“^ < p < 10“^, corresponding to 10“® < Ap < 10°. 
Most of our results are for N = 1, 000, but we also show 
data using N = 10, 000 to probe larger length scales. 

The total elastic modulus (bulk, shear), X{= K,G), is 
obtained as a sum of the affine, Xa, and non-affine, Xf^, 


components: X = X^ — X^ 12442811 . While Xa can 
be thought of as the value predicted assuming particles fol¬ 
low affine trajectories under an imposed deformation field. 
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quantifies deviations from this due to non-affine relax¬ 
ations. Yet, obtaining elastic modulus information in frag¬ 
ile systems can be problematic, especially when applying 
explicit deformation procedures. Here, we implemented 
protocols developed within linear response theory 
which avoid explicit deformation practices thereby allow¬ 
ing us to probe extremely close to the jamming transition. 

Two protocols were employed that essentially sample 
the vibrational normal modes of system: (i) The zero- 
temperature (T = 0) protocol (restricted to = 1,000) 
is formulated directly in terms of the dynamical ma¬ 


trix 112412511 . (ii) The finite-temperature (T > 0) protocol 
(for both N = 1000 and N = 10,000), samples mode 
vibrations by switching on a small temperature (T = 10“® 
to 10“^°) and thermally agitating the system At 

these temperatures and p > 10“^, particle displacements 
are 10“^ to 10“"^ [cr], and both protocols return consis¬ 
tent values. Technical details of numerical procedure and 
formulation can be found in Supplemental Material ll2^ . 
Here we highlight an important aspect of these protocols. 
Both procedures are accessible through current experimen¬ 
tal technologies at the colloidal and granular scales. In par¬ 
ticular, advances in particle tracking and resolution allow 
precision measurements of particle positions, used by co- 
variance matrix analyses methods 1134-^ . and the photo¬ 
elastic technique for particle forces I. 

To extract local information, the simulation box was 
divided into small subvolumes of size x Wy x w^, 
i.e. coarse-graining (CG) domains. In each CG domain m, 
we computed the local modulus, X'^ = AT™, G™, decom¬ 
posed into their affine (4) and non-affine {N) components. 
We then calculated the probability distribution function 
P{X™‘), from which the average X and standard deviation 
5X were obtained 1^ . 5X = 5X{p,Wx-,Wy,Wz) de¬ 
pends on both p and the size of CG domain, and quantifies 
the extent of fluctuations, whereas X = X{p) corresponds 
to the global value, independent of Wa (a = x,y,z) JUI]. 

Figure [T] shows the dependence on pressure, p, of the 
moduli and their corresponding fluctuations. The global 
X (p) are shown in the top panels. Fig. (Ua), (b), indicating 
that our technique is consistent with previous studies on 
similar systems |S[I1 that imposed explicit deformations. 
Since the pressure scales as p V ~ (~ 

for a = 2.5, Hertzian contacts), the scaling laws for X 
normalized by the effective spring constant ~ V" ~ 
ll^ . Xjkgff, are consistent with: 


iT/fceff ~ G/Ks ~ A(j)' 


0.5 


( 1 ) 


The middle panels. Fig. [He), (d), show the absolute fluctu¬ 
ations, 5X{p, Wx,Wy,Wz), where the CG domain is cubic 
of linear size, m ~ 3, and from which we And, 


5K/k,s ~ A(/.°, 5G /feeff ~ A(/.' 


0.27 


(2) 


More importantly, the bottom panels, (e) and (f), present 
the fluctuation data on a relative scale, bXjX, which gives 
the appropriate measure of the degree of heterogeneity. As 
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FIG. 1. (Color online) Elastic modulus dependence on packing 
pressure, p. The average (global) total, affine (A), and non-affine 
(A) values and corresponding standard deviations 

5X of the probability distribution P{X^), for A™ = A™ (left 
panels) and G™ (right panels). The CG domain is cubic of linear 
size w ~ 3. Lines are power-law scalings with p. The presented 
data were obtained using the T — 0 protocol with N = 1, 000 
and averaging over 10 different realizations at each value of p. 


A(j) —> 0 (p —> 0), 5K/K approaches a constant value, 
whereas relative fluctuations in the shear modulus grow as 

5G/G ^ Acj)-''^, zzg - 0.5 - 0.27 = 0.23. (3) 

We remark on two additional key features of Fig. [T] 
Firstly, for the bulk modulus the affine and non-affine com¬ 
ponents are quite distinct, such that the total bulk modulus 
is largely determined by the affine part only. Secondly, and 
in contrast to the above, the shear modulus components re¬ 
main close in value, so the scaling for total shear modulus 
is controlled by the gradual cancellation of affine and non- 
affine confributions. 

We now turn to a more explicit view of the spatial dis¬ 
tributions of iC™ and G™. Figure |2] presents the proba¬ 
bility distributions P{K'^) in (a) and P{G'^) in (b). We 
find that all the P{X™‘) are well-characterized as Gaussian 
over the entire pressure range, even down to the jamming 
point lH^]. But notice that although all the > 0, G"* 
can contain negative values. The fraction of these nega¬ 
tive shear modulus zones, P(G™)dG'", is 
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FIG. 2. (Color online) Probability distributions of (a) and (b) 
G™, and their relative fluctuations, (c) - K)/K and 

(d) G™ = (G™ — G)/G, for the range of p indicated in the legend 
of panel (a). Spatial correlation functions Gj^„ (r) (defined in 
main text) for (e) and (f) G™. The inset to (b) shows the 
fraction of negative G™ regions, Fn = f^m^Q P(G"^)dG"^, as a 
function of p. A close-up of P{K”^) is shown in the inset to (c). 
In (a)-(d), solid lines indicate Gaussians. In (e),(f), vertical lines 
indicate the CG length, r = Wa = in ~ 3. Data were obtained 
using the r = 0 protocol for = 1, 000. 


shown in the inset to Fig. Il^h). grows as p —)• 0, sug¬ 
gesting that there is a 1 : 1 ratio of stable and unstable 
regions |41| as the system becomes fragile ll3^ . Note the 
fact that our data appear to level off at the lowest pressure 
is likely a system size effect In Fig.|2c), (d), we plot 
P{K'^) and P[G^) of the fluctuations relative to global 
value, — X^jX. P{G^) broadens signifi¬ 

cantly as p decreases, which is quantitatively demonstrated 


by 5G/G in Fig. [Tlf) lEH], whereas variations in P{K'^) 
are rather small and insensitive to p, consistent with 5K/K 
in Fig. (He). 

In an effort to directly detect a correlation length as¬ 
sociated with these fluctuations, the bottom panels of 
Fig. llJe), (f) show the fluctuation spatial correlation func¬ 


tion, G^m{r) = ^X'"(r)X'"(0)y^i:”^(0)X’”(0)J), 

where we explicitly represent X™ as a function of posi¬ 
tion r, and () denotes a spatial average. Both the (7^m(r) 


(a) = 3 (b) dn, =2 



FIG. 3. (Color online) Dependence of SG/G on the CG length 
w for (a) dw = S and (b) dw = 2, as discussed in main text. 
Same symbols used as key in Fig.|2ta). Closed symbols are data 
using Wa = Ml ~ 3, i.e. same data as shown in Fig.[TJf). Lines 
are power-law scalings, (a) dGjG ~ and (b) ~ , 

consistent with dGjG ~ w~'^'“l‘^. Insets: same plots on log-log 
scales. For w < 6 data were obtained by the T = 0 protocol with 
N = 1, 000, and for w > 6 the T > 0 protocol with N — 10, 000. 


decay with the CG length r = m ~ 3 H, indicating 
that and G"* fluctuate randomly in space without any 
apparent correlation, which persists all the way down to 
the transition point. Thermal glasses HmiiSSl and 
disordered crystals ll^ Wh similarly exhibit random dis¬ 
tributions in their local moduli that are Gaussian. 

An alternative view to determining a possible charac¬ 
teristic length is through the dependence of fluctuations, 
SXjX, on the size of CG domain, x Wy x w^- We 
considered three different ways to change the CG domain: 
Vary, (i) Wx,Wy, equally, so that = Wy = = w, 

(ii) , Wy as Wx = Wy = w, keeping fixed w^ — 3, (iii) 
only Wx as w^ = w, keeping fixed Wy = w^ — 3. In (i), 
the CG domain is always cubic, whereas it becomes rectan¬ 
gular parallelepiped in (ii), (iii). We define the dimension 
dw of CG domain; dw = 3, 2, 1 for (i), (ii), (iii). As 
we have seen so far, X"^ is a random variable, following 
a Gaussian P{X'^). Thus, within the framework of a sum 
of random variables ll^ . we obtain the scaling law with 
respect to CG length w: 


6X/X ~ 


(4) 


Figure [3] shows the m-dependence of 6GjG at several dif¬ 
ferent p, for = 3 in (a) and dw = ‘2 in (b) (see for 
dw = 1). For all pressures, 6GjG ~ w~^'‘^, ~ w~°'^^, 
~ y;-o.47 Pqj. _ 2 ^ 2^ 1, respectively, which all confirm 
Eq. (lU). We obtained the same result in 6K/K. The same 
power-law dependence on w has been reported for glasses, 
with exponent 1 in =2 lld^ and 1.5 in = 3 lE^ . 

Combining the scaling results for bGjG (Eqs. (l3]l 
and dU)), expresses that relative fluctuations in shear modu¬ 
lus are suppressed over sufficiently large w. This supports 
the existence of a characteristic length, above which 
fluctuations become negligible. Specifically, we define 
as w at which we see a fixed value, ccq, of bGjG for all p 
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(c) p = 4 X 10 w ~ 6.5 



X 


FIG. 4. (Color online) Spatial maps of local shear modulus fluctuations, G™ = (G™ — G)/G, within a fixed x-y layer, for dw = 3. 
(a) Large pressure p = 4 x 10“^ and small CG length w ~ 3, (b) small p = 4 x 10“^ and small w ~ 3, and (c) small p = 4 x 10“^ 
and large w ~ 6.5. Data were obtained using the T > 0 protocol with N = 10, 000. Additional snapshots shown in Supplemental 
Material 1^ . 


or Ai/), i.e. we determine as (5G/G = aQ{w/^G) 
which gives ll^l^ 

u^ = UGl{d^l2). (5) 

The idea of the length ^g associated with growing 6GjG 
is best visualized in Fig. |4l which shows the local fluctu¬ 
ations of shear modulus (for = 3) as follows: Pan¬ 
els (a) and (h) of Fig. |4] compare modulus maps of G"* = 
(G™ — G)/G for a slice through two packings at two dif¬ 
ferent p, at the same w ~ 3. In relation to Fig. [Sja) 
(dw = 3), these two points lie at different values of 6GjG 
along a vertical line at m ~ 3, that intersect the respec¬ 
tive p curves. At this value of w, the two systems appear 
very different. Far from (pc, Fig- Ha) (p = 4 x 10“^), 
the system appears quite uniform, and fluctuations are sup¬ 
pressed. Whereas, close to cpc. Fig. |2h) (p = 4 X 10“®), 
we observe large-scale, spatial fluctuations. For the system 
closer to pc (small p), fluctuations become suppressed at 
the larger w = 6.5 (Fig. Uc)), so that the map resembles 
more compressed system at the smaller value of w. This 
corresponds to drawing a horizontal line across Fig. [Sja) 
at the same value of bGjG connecting the two curves at 
different p. 

In conclusion, we found that the differeces between bulk 
and shear moduli fluctuations, as the jamming point is ap¬ 
proached, are caused by the non-affine components. Rel¬ 
ative fluctuations in the bulk modulus become insensitive 
to packing pressure as !\p 0. Whereas, shear modu¬ 
lus fluctuations increase as, bGjG ~ which leads 

to the identification of a lengthscale, ^g ~ A(/>“‘^^. For 
CG dimension, = 3,^ ~ 0.16, a value distinct from 
any previous study ifToU^ Isilj^ . ^g corresponds 
to a scale above which the elastic properties coincide with 
those of the bulk system, while below, the local mechan¬ 
ical properties deviate from macroscopic behavior. It has 
been proposed that a continuum elastic description breaks 
down below a scale, £c ~ in two dimensions ll2flll . 

consistent with our ^g for = 2, and can derive from the 
transverse component of elastic response J^, which are 
controlled by shear modulus fluctuations. 


At the same time, however, we also found that the lo¬ 
cal elastic moduli randomly fluctuate without any appar¬ 
ent correlations. This feature seems to be general for a 
wide class of disordered materials, thus further promot¬ 
ing the idea that granular-like particle systems present a 
model state for examining mechanical properties of disor¬ 
dered materials. Curiously, the randomness in local moduli 
persists down to the transition point and is different from 
the distribution of contact forces, which becomes more ex¬ 
ponential closer to pc 10] and is therefore more sugges¬ 
tive of spatial correlations. Such random fluctuations in 
the moduli may come from the coarse-graining procedure 
and/or the random distribution of particle contacts, which 
is a topic for future investigation. 
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